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NOTATION 


A, B, C, 
a, b, c, d, 


a 



k 0 

K 

L 

m 

n 

q 


= principal moments of inertia 

= numerical coefficients, Eq. (64) 

= correctional factor, Eq. (27) 

= drag coefficient 

= lift coefficient 

= pitching moment coefficient 

C , C = stability derivatives 
m ^ m a 
“ q 

= coefficient, Eq. (48) 

= coefficient, Eq. (54) 

- acceleration due to gravity 

= correctional coefficient, Eq. (41) 

= radius of gyration in pitch 

= A - C 
B 

= coupling coefficient, Eq. (54) 

= mean chord, characteristic length 
= mass of vehicle 
= angle-of-attack frequency 

= angular velocity in pitch relative to the earth 


L 

2u 0 q 


= radial distance from center of earth 
r 


s 

S 

t 

ft 

ft 

T 

u 0 


= — jz °. .. , speed ratio 

^ go^o 

= reference area 
= time 



V 


u 0 

= thrust, period 
= reference circular speed 


v 



w 

V 
X 
x 
a 

P 

V 
6 

e 

■n 

e 

x 

H- 

e 

p 

cr 

T 

CO 

Subscript 
( )o 
< >i 


= vehicle weight 
= speed along the orbit 
= vector state, Eq. (13) 

= cos T 

= angle -of-attack 
= density coefficient Eq. ( 11) 
= flight path angle 



= small perturbation, eccentricity of orbit 
_ 

L 

= angle of pitch 

= characteristic value 

PoSu 0 
2m g 0 

= density coefficient Eq. (24) 

= air mass density 

= non-dimensional density gradients Eq. (10) 

= - w -^t 

u 0 

= central range angle (Fig. 1), function, Eq. (66) 
= phugoid frequency 


= reference flight path 
= perturbed quantity 
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ABSTRACT 


This paper presents an analytical study of the longitudinal dynam- 
ics of a thrusting, lifting, orbital vehicle in a nearly circular orbit. 

The translational motion is composed of a non-linear oscillation, or 
phugoid, and a spiral mode which results in either decay or dilatation of 
the orbit depending on the perturbed initial conditions. The non-linear 
effects on the phugoid period and damping are small in the altitude range 
considered. Elements of the orbit such as radial distance, velocity, 
and flight path angle were obtained explicitly as functions of time. The 
behavior of the variations of these elements is correctly predicted. 
Explicit expressions for period and damping of the angle-of-attack mode 
were derived. It is shown that a critical altitude may exist at which the 
phugoid mode and the angle-of-attack mode have nearly equal periods. 
Near this resonance altitude linearized solutions are no longer valid 
and a study of the non-linear equations shows that there is a strong in- 
teraction between the translational and the rotational modes resulting in 
a switching of the two frequencies of oscillations. 
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NON -LINEAR LONGITUDINAL DYNAMICS 
OF All ORBITAL LIFTING VEHICLE 

INTRODUCTION 

In Ref. I, Etkin has presented a very enlightening study of the small-pertur- 
bation dynamics of a satellite vehicle in a nearly circular orbit. He linear- 
ized the equations of motion and solved the resulting fifth-order system 
numerically. It was found that the linear solutions contain two oscillating 
components which can be identified with the classical phugoid and short-period 
modes and a new spiral mode. Etkin then showed by direct numerical calcu- 
lations that, with hypersonic speeds at a flight altitude where the gravity 
torque predominates over the aerodynamic torque, the so-called short-period 
oscillations can develop a period that is longer than that of the corresponding 
phugoid. At the critical resonance altitude where the two periods are nearly 
equal, Etkin f s linearlized solutions are no longer valid because the ampli- 
tudes of the oscillations are large. 

More recently, E. V. Laitone and Y. S. Chou made a theoretical analysis of 
the same problem. 2 Their analytical solution of the linear equations are in 
excellent agreement with Etkin f s numerical calculations. 

In this paper, we extend Laitone and Chou' s investigation to include non- 
linear effects in the longitudinal dynamics of the orbiting vehicle. The equa- 
tions of motion of a thrusting, lifting vehicle in a nearly-circular orbit are 
integrated directly in matrix form, using a perturbation technique. It will be 
shown that, for the phugoid, or trajectory mode, Etkin f s new spiral mode is 
the familiar secular perturbation of a vehicle flying in a resisting medium. 

The proper phugoid motion, with damping, is described by oscillatory terms 
with diminishing amplitudes. Through explicit formulas, the asymptotic be- 
haviors of the variations of the elements of the orbit are correctly predicted. 
Furthermore, a very simple formula yields a value for the altitude where 
variations in velocity change sign. This value for the velocity inversion alti- 
tude is accurate to within 10 feet of the numerically computed value. For the 
short-period mode, which we shall refer to as the angle-of-attack mode, ex- 
plicit formulas are derived which yield accurate values of the period and 
damping at all altitudes. The resonance altitude where the two oscillatory 



modes have equal periods is obtained by solving a very simple equation. 
This value for critical altitude is also accurate to within 10 feet of the exact 
numerical solution. 

THE EQUATIONS OF MOTION 



Fig. 1 Axes System and Nomenclature 


If we use an axes system that is always tangent to the flight path as shown in 
Fig. 1, the motion of a lifting vehicle with constant thrust is governed by the 
system of equations* 


dV _ 

T 

— cos a - 
m 

P SC D V2 

dt 

2m g 
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dt 
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0 - y + a 


*A11 symbols are defined in notation section. 
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The first two equations are, respectively, the drag and lift equations along 
the tangent and normal to the flight path. The constant thrust, T, is set for 
a reference circular orbit of radius r 0 with constant lift and drag coefficients. 
Furthermore, the thrust is small enough that the mass of the vehicle can be 
assumed constant. Hence, along the reference orbit, where y = a = 0 

T=i P 0 SC D( U 0 2 


1 


up 

go r o 


PoSC 

J-'O 
2m g 0 


2 

u 0 


( 2 ) 


C 

m 0 


0 


qo - 


uo 

r 0 


The first term on the right-hand-side of the pitching moment equation, Eq. 
(1-3), expresses the restoring aerodynamic torque, while the second term 
corresponds to the gravity torque. The last three of Eqs. (1) are kinematic 
relationships. The mass density, p, of the atmosphere is solely altitude 
dependent. For computational purposes the atmospheric data used were ob- 
tained from a polynomial representation of the 1962 U.S. standard atmos- 
phere as presented in the 1966 U.S. standard atmosphere supplements. The 
lift and drag coefficients, C T (o') and C (a), are functions of the angle-of- 
attack only, while the pitching moment coefficient q) depends on both 

the angle-of-attack and the angular velocity in pitch relative to the earth. 

To write the equations in non-dimensional forms, let 


V(t) = u 0 u(t), r(t) = r 0 r(t),q(t) = ^ q(t) 
p(r) = p 0 p(r), g(r) = go(^r) = go ( ^ ) 


Furthermore, 


to the first order 
C D («) 

C L<“> 



+ C a 


+ C T a 

LlQ, 


(4) 
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c (a, q) - C + C a 

m xxi q xxx qj 


L 

2u 0 


C (q -q 0 ) 


m 


where C_ , C T and C are the drag, lift and pitching moment coefficients 
Dq -Lq m 0 

along the reference flight path. 


With the definitions 


t = &t. k 2 = 

) ■(•) 
d t 


U 0 y m 


_ PoSuo _ 

11 ' 2mg 0 ' k » - 
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2r 0 

L 

2 

2 - U 0 

, _ ( L \ 2 


S “ , 

go r o 
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(5) 


the non-dimensional equations of motion can be written as 

A A A 2 ^ \ 1 

u = |J,C cos cx - ppu C (tt) - — sin Y 


D 0 


A. . AA 2 . / 1 S* U* 

UY = ^ C D 0 Sma ^ pU C L a 


q = Zufipu'C^or, q) 


3k 

A, 

2t| r 3 


^ sin 20 


cos Y 


( 6 ) 


A 

* ? A 2 U 

0 = s rjq + s — cos y 


A 2A 

r = s u sin y 


0 - v + 


y -r a 


We see now that, for prescribed initial conditions, integrating the system 
of Eqs. (6) is a formidable task. The task is eased, however, by decoupling 
the equations. To this end, we shall assume that the motion of the vehicle 
around its center of mass has negligible effect on the orbit, that is, on the 
trajectory mode. This assumption is justified since the vehicle is small in 
comparison with the dimension of the orbit and since the aerodynamic forces 
in the altitude range considered are also small. This statement of the so- 
called limited problem has been used in celestial mechanics to study the 
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libration of the moon. Based on this decoupling, we shall first calculate the 
trajectory mode, and then use the results to integrate the equations govern- 
ing the angle -of- attack oscillations. The coupling effects will be examined 
in the last part of the paper. 

PHUGOID OSCILLATIONS 

The phugoid, or long period oscillations, occur at nearly constant angle-of- 
attack. Hence, by taking a ~ 0, we have the non-dimensional equations 
which govern the phugoid mode. 

u = fiC- (1 - pu ) - — sin v 
D 0 a* 

A ■ /n Ex AA 2 / 1 s 2 U* \ ^ 

uy = (1 - s )pu - ( — — J cos v ( 7 ) 

\ r 2 r ^ 

A 2A 

r = s u sin y 

We note that s is the ratio of the circular velocity to the orbital circular velo- 
city without drag. At very high altitude s -** 1 and p -► 0 and the first and the 
third of Eqs. (7) reduce to 

A 

u = - 

a a 
r = u 

From these, we have the energy 

i A ? 

i u2 - 

In general, the force field is not 
by atmospheric drag. 

If we allow only small departures from the equilibrium flight path, we can 
express the orbital parameters in terms of small perturbations 


A sinY 

r c 

sin y 

integral 

1 

— = constant 
£ 

conservative since the energy is dissipated 
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G( £) = 1 + uj ( £) 
r(£) = 1 +r l( 1) 
Y(*> = Yi(*> 


( 8 ) 


Also, by expanding the mass density p(r) in a Taylor's series in r t , we have 


where 


A/A. . A A2 

p(r) = 1 + o-, rj + o- 2 rj + 


= /dp \ £o „ = i/d 2 p\ r| 

1 \dr J 0 p 0 ' 2 ^\dr 2 J> p 0 ' 


(9) 


(10) 


Further, let 


u 2 = (1 - s 2 )(-cr 1 s 2 + 2) + s 4 
P = -is 2 [(l -s 2 )(a 2 - 1) -2] 
t “at £ ( ) = ( )' 


( 11 ) 


To the second order of smallness in the perturbations, we can rewrite the 
system of Eqs. (7) in matrix form 


X'(t) = AX( t) + B(X) 


where 


(12) 


X(t) 


“A 


' 2|,C D 0 - 1 

Ul 




1 

Z - co 2 

Yi 

, A = - 

CO 

2 0 2 — 

s 2 

A 





0 s 2 0 


(13) 


B(X) = - 
w 


P C D 0 ( ^f + 2(ri ^ 1 ^ 1 +cr 2 ^i) + 2 Y 1 ^ 1 
-Uj +-|(1- s 2 )yi - |r - -4 (w 2 - 2 + 4s 2 )r , 1 


s 4 Uj Vl 
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Let € be an initial perturbation of any one of the variables, say r x . We can 


then assume series solutions of the form 

A A , A 

u i = €u n + e u iz + * • ■ 
Yi = € Yn + e z Yiz + • • • 

(14) 

ri = €r n +€ 2 r u +• • • 


or 

X = €Xj + € 2 X 2 + • • • 

(15) 

where the definitions of the column vectors Xx , X 2 , . . 

. are clear from Eq. 


( 14) . By substituting into Eq. (12) and equating coefficients of like powers 
in e we have the system of differential equations 

Xi = AXj ( 1 6) 

X' 2 = AX 2 + B(X t ) (17) 

Since A is a constant matrix we can immediately integrate Eq. (16) to have 


__ / . At__o 
X l (t) ”6 X; 


( 18 ) 


where 


X? = Xj(0) = e -X(0) 


(19) 


With this value of Xx(t) we can express B(X x ) as a function of t and consider 
Eq. (17) as a linear system of differential equations with a forcing term. 

By integrating we have 

T 

X 2 (t) = e AT X® + ye A(T_t) B(t) dt (20) 

0 

But 

X® = X 2 (0) = 0 

Hence 

T 

X 2 (t) = j* e A(T_t) B(t)dt (21) 

0 

and the solution to the second order, of the phugoid mode, is 
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, _ At „0 , _2 

x(t) = ee Xi + e 




-t) 


B(t)dt 


( 22 ) 


Linear Solution 


The integration of Eq. (16) gives the linear solution. The system has the 
characteristic equation 

2^CL _ 2pC. 


X 3 + 5o\ 2 + X - 


^(q = o 

\w J 


where 


£ 2 = s 2 [-(l+cri)s 2 +2] - -w 2 - cr i s 2 +2 


(23) 


(24) 


In general, the characteristic equation has a pair of complex conjugate roots 
corresponding to the phugoid oscillations and a real root corresponding to 
the spiral mode. The last term of the characteristic equation induces the 
spiral mode and is a small quantity. Let 


2pC 


Do 


(t 


spiral co 

where a is a quantity to be determined. By substituting into Eq. (23) we 
have 

\2 


(25) 


Using Lagrange f s expansion 3 we have for the value of a 

\ CO / \C0 J \ \(0 J J \ co J \co / \ \ CO / J\ \co 

- (— ^^-^ 6 (l+( / -^(l + 2r- N ) V5 + 6^ 

\ CO J / \ V 00 J / V J J\ \W 


( 26 ) 


a = 1 


(27) 


Then by factorizing the cubic equations, Eq. (23), we have for the phugoid 

mode ^ 

2pC 


X 2 + 


a jX + a - 0 


This gives 


Real (X , . ,) 

phugOLd 


/ /£ V 

[l+[-) 

\ \co / 

+ (t 

CO \ \co 


(28) 
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Im (X. ) 

phugoid 




? \Ti 


( 29 ) 


With the roots calculated we have for the spiral mode in real time 

0.69m /~>?r / 2 M-C, 


double p 0 SC u 0 
JJo 


(f)K^Xt>(K§>)] - 


For the phugoid mode the damping is given by 

1.38m 1 


tolf ' « sc d, u ° T^TJ, 


In real time, the phugoid period is 

2vu- r /^n-v 


T = 


ugo 




--2 


(3 1) 


(32) 


In the last equation, by taking the bracket equal to unity we have a result that 
is identical to Laitone and Chou's Eq. (1. 5) 2 . At very high altitudes, co tends 
to unity and p tends to zero. The phugoid period asymptotically tends to the 
circular orbital period. In reality, after a perturbation has been applied in 
a vacuum, the vehicle will go into a slightly elliptical orbit. Thus, T should 
tend to this elliptic orbital period. This correct orbital period appears only 
when we consider non-linear terms. 


In general, the quantity 


(^Xi 


\ 2 


) in the expansion of a, Eq. (27), is 

/ 


small and setting a = 1 gives a very good approximation. This expansion 
gives the roots of the characteristic equation explicitly to the desired de- 
gree of accuracy. In our derivation, the damping term for the phugoid is 


exp 


p 0 u 0 SC 
E_o 

2m 



(33) 


where a is explicitly given by the series expansion, Eq. (27). Using a = 1 
we have Laitone and Chou # s Eq. (3.9) for phugoid damping. Hence, besides 
the additional spiral mode obtained, the above results greatly improve the 
already accurate formulas for phugoid oscillations derived in Ref. 2. 
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Now let 


Vugoid = ± iwi 

= ^2 


(34) 


’spiral 

With the roots calculated, we have, to the first order, for the elements of the 
flight path 

r = 1 + €Cie^ 2T + ee^ lT [C 2 COsoj!T + C 3 sinoji t] 

Y = ir' C^e + e -—2 e 1 [(W1C3 +\iC2)coscoii (wiC 2 ~ X1C3 )sintoiT] 


u = 1 + T~r [(u 2 - Z) + gj 2 \ 2 ]c 1 e 


(35) 


€ \l 

+ T7 e 
2s 


T {[^ (0j2 ~ ^ +c ° 2 1 - <j^i ) ]C 2 + 2w 2 w 1 \ 1 C 3 J cos co 1 t 

[(to 2 - 2 ) +cj 2 (K j - coi ) ]C 3 ~ 2 oj 2 coi^iC 2 J s in 00 ^ t "j- 


where the Ch are constants of integration. 

From these expressions we have the following interesting remarks: 

1. In each of the variations of elements such as radial distance, flight path 
angle, and velocity, there are two components. One component is oscilla- 
tory with diminishing amplitude and it tends to circularize the flight path. 
The other component is aperiodic and divergent. This component is due to 
the offset effect between the thrust and the drag and induces a secular varia- 
tion of the elements of the orbit. 

2. For the vehicle considered in Ref. 1, coi > 1 above about 140, 000 ft, and 
o>i -► 1 as r 0 -* 00 , thus the effect of drag is to shorten the phugoid period. 

3. In the expression for r, since \ 2 > 0, the divergent mode tends to de- 
crease the radial distance if the initial perturbations are such that €Cj < 0. 
On the contrary if €Cx > 0 the radial distance will increase with time. 
Furthermore, under the constant thrust application, with decreasing drag, 
the vehicle will move outward following a spiral. 

4. From the expression for the flight path angle we can see that it varies in 
the same direction as the radial distance. 
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5. On the contrary, the velocity varies as the radial distance if and only if 


> r fa (36> 

If the inequality is satisfied the velocity will increase as the vehicle is 
spiraling out and decrease if the vehicle is spiraling in. The inverse is true 
if the inequality reverses. To calculate the altitude where the velocity in- 
version occurs, since is small, we can use the equation 

w 2 = 2 (37) 

The Eqs. (3 5) are derived in terms of the non-dimensional time t, and the 
frequency of oscillations u>i is given by Eq. (29). In terms of the non-dimen- 
sional time, the frequency of oscillations is o^to which can be approximated 
by co since co! is near unity. The square of the linear phugoid frequency, <o z , 
is plotted versus the altitude in Fig. 2. 

We see that to 2 is large at low altitudes and tends asymptotically to 1 when 
the altitude increases indefinitely. When co 2 = 2 the velocity inversion occurs. 
More explicitly, using the definition in Eqs. (11) for co 2 , we have 

(W/S) g r 0 

~C = ‘ ~ 2 -< 2 Po + P f o r o) (38) 

Lo 

where subscript, s, denotes the condition at sea level and p f 0 is the density 
gradient evaluated at r 0 . The left-hand side of the formula above is a char- 
acteristic of the vehicle and the right-hand side is solely dependent on the 
characteristic of the atmosphere. Fig. 3 is a plot of Eq. (38) as a function 
of the altitude and the graph can be used to. determine the altitude where the 
velocity inversion occurs for any given vehicle. For example, the vehicle 

considered in Ref. 1 has (W/S) /C T = 600 lb/ ft. 2 Thus, the critical alti- 

s Ljo 

tude for velocity inversion is 321, 000 ft. 

Eq . (38) gives the critical altitude to within 10 ft. of the value computed 
numerically from the exact linear equations. 
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Non-Linear Solution 


We have seen that the next higher order component of the vector state is 
e z X 2 . If B(t) SB1X,), then 

T 

X 2 = Je A(T_t) B(t)dt (39) 

0 


and, A and B are given by Eq. (13) and X 1 represents the linear solution 
found in the preceding section. Hence the computation of X 2 (t) is straight- 
forward. Here we need only to have an idea of the order of magnitude of the 
second order term. By the form of the vector B(Xi), we can see that, if the 
perturbation is small, the contribution of € 2 X 2 is negligible. This has been 
shown in a numerical study by Rangi. 4 But in the expansion, Eq. (9), of the 
mass density of the atmosphere, cr 2 is of the order 1 0 5 . Hence, if e>0{10 ) 

we should include the second order gradient effect of the air mass density <x 2 
which appears explicitly in B(Xi). By the same cons iderat ion, terms which 
need to be retained in B are 




D 0 A2 
— (f z r n 


B(X x ) = 


co 


J-JL £ 2 
rn 


(40) 


Hence, the non-linear numerical analysis of Ref. 4 is not valid for large per- 
turbations, since the author has neglected all second and higher order terms 
of the air mass density. For large perturbations, as were considered in 
Ref. 4, the prime contributing non-linear factor in the phugoid and spiral 
trajectory is the variation of the mass density of the atmosphere. The per- 
turbed trajectories for these cases are highly eccentric orbits and it is not 
correct to assume a linear variation for the air mass density. (See Appendix C) 

The value of o- 2 considered above is somewhat too large because of an inverse 
polynomial representation of the atmospheric mass density. This resulted 
from a curve fitting analysis which can give a wrong value for a truncated 
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series at a certain altitude. In trajectory analysis a better approximation is 
usually found using an exponential atmosphere. However, this atmosphere is 
not suitable for a dynamics stability analysis in which a series expansion of 
the air mass density is required. 

In this paper, we just want to call attention to the effect of the second order 
atmospheric gradient. For an accurate second order analysis, the coeffi- 
cients cr i and c r 2 in the "parabolic representation" of p, Eq. 9, should be 
averaged for each limited altitude range considered. 

Asymptoti c Beh avio r of P h ugoid Period 

To the first order we have seen that the phugoid period tends to the circular 
orbital period when the altitude of flight increases indefinitely. In reality, 
the phugoid period tends to the perturbed elliptical period. To find this -cor- 
rect orbital behavior we use a new time variable, t, such that 

t = t( 1 + hj€ + h 2 € 2 + * * *) 

where h x , h 2 are constants to be determined. By substituting into Eq. (39), 
neglecting drag terms, and requiring a periodic solution for X 2 , we can 
easily find that 

_ (w 2 - 2) 2 - 2(1- s 2 )(w 2 -2) +4s 2 <(3- 2) 

- hl (41) 

Hence, to the first order, an asymptotic expression for the phugoid period 
is 

T = (1 + eh! + ■ • •) (42) 

wg 0 


ANGLE OF ATTACK OSCILLATIONS 


The angle-of-attack mode is governed by the system 


dq 

dt 


pV 2 SLC m (*, q) 
2B 


3g( r) ( A - C) 
2r B 


sin 20 


d0 V 

_ = q + - cosy 


0 = Y + a 


(43) 


Ik 



The elimination of 0 and q results in a non-linear equation in a 

d 2 o; , T cos a da, pSV dC ]J a) 3g (A - C) 

312 " + rr~ 337 + 3 — — 3i — + o 5 — — cos 2 Y sin 2a 

dt 2 mV dt 2m dt 2r B 

+ ■ sin 2 y cos 2a + -^- •— (C T {a) cos a + C^a) sin a) 


m 2m 


X g C T si 

+ — ^t( 2 sinv sina; + cos v cos a) - ( — ) — 
mV 2 1 \ m ) t 

- (H ) ! c d 


in 2a 


pSC (a) 

,(a)C L (a)V 2 - — 2 ^— S cos y 


- £ ® !ic n, ( “' q) + l? sln2 '' + is- af -^ sin21 ' = 0 


(44) 


Linear Solution 

A zero-order solution can be obtained easily by considering small oscilla- 
tions of the angle-of-attack along a circular orbit. Assume 


or « 1. Y = 0, = 0, V = u 0 , r = r 0 


(45) 


Then Eq. (44) reduces to the linear equation with the non-dimensional time, 
t, as the independent variable 


+H (C +C - 26C )«' +T-(3k 0 -2|xri6C )a = o 

to Do m A w z ^ m a 


(46) 


1 A 

q 

To the order ^ we have for the angle-of-attack mode 

oj co 


(47) 


where 


“■itC^+C, 26C ) 


'Ni 


'Do 


a 


m A 

q 


(48) 


n* = s ( 3k 0 - 2 J 1 T 16 C ) 

m ff 


The solution, Eq. (47), is identical to Laitone and Chou's result 2 and is in 
good agreement with the results of the numerical study in Ref. 1. 

Explicitly, in real time, the angle-of-attack period is 


15 


1 

“ ~7 


T = 


f^7 l 3ko ‘ 2 ^ 6C 


m 


a -* 


and the damping is given by 


M-go 


exp I - ““ (C: + C T - 26C )t 

K 1 2u 0 D 0 L m A J 

“ q 


(49) 


(50) 


Resonance Altitude 

For practical values of vehicle parameters, there is an altitude where the 
two frequencies, to and n, can be equal. This altitude, called the resonance 
altitude, is found by solving the equation to = n, and for simplicity, setting 
s 2 = 1. 


Explicitly, we have 


2k 2 (3k 0 - 1)(W/S) 

_X £ 

LC 


m 

a 


Po^ogg 


(51) 


where subscript s denotes the condition at sea level. The left-hand side of 
the formula above is a characteristic of the vehicle and the right-hand side 
is solely dependent on the characteristic of the atmosphere. Fig. 4 pre- 
sents the variation for the earth's atmosphere as a function of the altitude 
and the graph can be used to determine the resonance altitude for any given 
vehicle. For example, for the vehicle considered in Ref. 1, the characteristic 
value on the left-hand side of Eq. (51) is 5.2281 x 10 4 lb / ft 3 and the resonance 
altitude is therefore 49 2, 300 ft. The value of the resonance altitude com- 
puted from Eq. (51) is correct to within 10 ft compared with the exact value 
from a numerical analysis using Etkin 1 s equations. 
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Fig. 4. Resonance Altitude 


Coupling Effects 

We have assumed that the coupling effects between the two oscillatory modes 
are negligible. This is, in general, true, below and above the resonance 
altitude. However, near the resonance altitude, strong coupling effects are 
evident and the linear solutions, Eqs. (29) and (47), are no longer valid. 

In his study 1 , Etkin found that the Phugoid period increased with altitude, 
tending asymptotically to the orbital period. He also found that, for his ve- 
hicle, the period of the angle-of -attack mode increased with altitude and 
crossed over the Phugoid period at about 490, 000 ft, tending asymptotically 
to infinity at about 505, 000 ft. This behavior is also predicted by the un- 
coupled oscillatory mode frequencies, Eqs. (29) and (47), and agrees with 
one’s physical intuition. However, a close inspection of the region near the 
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resonance altitude, using Etkin f s linear equations, shows that this descrip- 
tion of the behavior of the mode periods is in error and reveals, instead, 
the phenomenon shown in Fig. 5. There is, in fact, a "switching" of the 
modes of oscillation at the resonance altitude instead of a "crossing"! 

The uncoupled damping constants in Eqs. (29) and (47) show that they will 
remain negative at all altitudes. However, Fig. 6 shows that near the reso- 
nance altitude the Phugoid damping constant becomes positive while the 
angle -of -attack mode damping constant remains negative. Above the reso- 
nance altitude the "switched" angle-of-attack mode damping constant is 
positive while that of the "switched" Phugoid is negative. 


A factorization of the fifth-order characteristic equation, of the linearized 
coupled equation of motion, that takes into account the effects described- 
above, has been obtained by Dobrzelecki. 5 He found that for the damping 
constants 


Real (\ . ) = - 

phugOLd 


r m- c , 


CO 


■(1 + 


ft)' 


11 ■ IC N, *vl ISi) 


Real (\ 


attack 


M-C 


N, 


+ K-^(C 


4co N i 




(53) 


where 


K = 


2[s 2 (1+s 2 ) + 3k 0 s 2 - co 2 ] 


_ „ iL 

w 2 [l - (*) ] 


CO 


C n 2 “ 2 


= 4[C n + C T + 26C ] 

2 D 0 L m A 

“ q 


(54) 


and the coupling coefficient K is positive below and negative above the reso- 
nance altitude. 

Eqs. (52) and (53), together with the equations for co and n, correctly and 
accurately predict the values of the damping constants and frequencies of 
the oscillatory modes, to within a few thousand feet below and above the 
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Fig. 5. Frequencies Near Resonance Altitude 



ALTITUDE (FT XIO" 3 ) 


Fig. 6. Damping Constants Near Resonance Altitude 
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resonance altitude. Figs. 5 and 6 clearly show that near the resonance alti- 
tude one has to consider the non-linear equations of motion to predict cor- 
rectly the crossing of the periods. 

Non-Linear Angle-of-Attack Frequency 

At high altitudes, especially near the resonance altitude, the amplitude of 
oscillation of the angle-of-attack may be large. The mode is governed by 
the non-linear Eq. (44) . A simplification may be made by observing that the 
aerodynamic forces are small and, thus, we can consider the orbit as an ex- 
act Keplerian ellipse for the duration of one or several revolutions. Further 
more, by considering a circular orbit and neglecting quantities of the order 
p 2 and the negligible damping, we have an equation with r as the independent 
variable. 


M , 3s k 0 . 2ps t]6 _ 

a' + sin a cos a - z ■■ C 

CO 2 Cl) 


m 


= 0 


(55) 


This non-linear equation is similar to the equation derived by V. V. Beletski 
(Ref. 6, p. 72) with one difference. Here we have a constant low thrust 
while Beletskii explicitly assumed a drag-free orbital period for the trajec- 
tory. 

For a slender vehicle with a conical surface of attack, an approximate ex- 
pression for the pitching moment coefficient is 


C = C sin o' cos a ( 56) 

mm 

a 

where can be evaluated by using the simple Newtonian impact 

a a __ 

theory for moderate angle-of-attack. Accurate values for C can be ob- 

a 

tained from wind-tunnel measurements. With this assumption, the non- 
linear equation for the angle-of-attack is 

-2 

a " + — 5” sin a cos a = 0 (57) 

Gu 

where 

n 2 = s 2 (3k 0 - 2 |jji 6C ) « n 2 (58) 

a 
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Eq. (57) can be integrated easily using the theory of elliptic integrals or the 
so-called Lindstedt method for obtaining periodic solutions to non-linear 
equations (Ref. 7, p. 141). The period of oscillation for large angle-of- 
attack obtained, in real time, is 


T = |3>(3k 0 - 2^6C m )<1 + + ■ ■ •) (59) 

go a 

where <*0 is the initial perturbed angle of attack. We notice that resonance 
does not appear for circular orbits. 


Eccentricity Oscillations 

If the orbit is elliptical, the result is qualitatively different since the co- 
efficients of the non-linear equation are periodic quantities. The most sig- 
nificant effect is the forced oscillations due to the non-vanishing eccentricity 
of the orbit. This effect will give rise to a possible resonance. Let us con- 
sider a slightly perturbed Kepler ian ellipse about the reference circular 
orbit, then 


= r 0 (l - e 2 ) 

1 + € COS T 


(60) 


where € is small. To the order € we have 


r = r 0 (l -€ cos t), V 2 = g 0 r 0 (1 + Ze cos r) 


V = (g 0 r 0 ) 2 (l+e cos t), g = go (1 + 26 cos t) 

GO . GO ., _ . 

sin \ = — € sin t, cos y - - (1 - Ze cos t) 

S S 

P = pofl-Ccrj cost), T = - “ Q( ° t 

u 0 


(61) 


Neglecting both damping and quantities of order p. 2 , and taking go = s = 1, we 
have 

a" + [(3k 0 - Zpr|6C ) - £(3k 0 - 2 ( 0 -! - 2) \xr)6C ) cos Tjsin a cos a 
a a 


+ 3£k 0 sinT cos Za + {jlC cos a 

Do 

= lxC^ (1 - €01 cos t) - € sin t - u.C T ai€ sin t 
•Dq Do 


( 62 ) 
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To display the effect of eccentricity oscillations, let us consider the case of 
small angle-of-attack. This gives us a Mathieu equation with forcing terms 

a" + (a 2 - eh cos t )a = -€(c sin t + d cos r) ( 63) 


where 


a 2 = 3k 0 - Zp.rjSC^ 

a 

b = 3k 0 - Z(cri - 2) uni 6C 

m 

a 

c = 3k 0 + 1 + jaC r 


(64) 


d = ^ C D 0 ai 


The homogeneous Mathieu equations have a periodic solution which can be ob- 
tained with classical methods. Here we seek a particular solution for the 
non-homogeneous Eq. (63) with each of the two forcing functions. For ex- 
ample, consider 

a" +(a 2 -€b cos r)a = - Cc sin t (65) 

and assume a particular solution of the form 


- - 6c s in t<(>( t) 

By substituting into Eq. (65) we have an equation for <j> 

(1 - x 2 ) - 3x + (a 2 - 1 - €bx)(j) = 1 

dx dx T 

where 


x = cos T 


( 66 ) 


(67) 


( 68 ) 


A particular solution of Eq . (67) is sought as a series in terms of the small 
parameter C. Let 

<H x) = 4> 0 + e 4>! + e 2 ^ + . . . (69) 

Then by substituting into Eq. (67) and equating terms of the same order of 
magnitude we have 
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(1 - x 2 )cf >0 - 3x4>0 + (a 2 - 1)4> 0 = 1 
(1 - x 2 )<j>i + (a 2 -l)4>i = bx<() 0 


( 70 ) 


where the prime here denotes differentiation with respect to x = cos t. The 
system can be readily integrated to give 



b cos T 

^ 1 ( a 2 - l)(a 2 -4) 


(7 1) 


It is clear that we can do the same with the forcing term -€d cost, and we 
have a particular solution for the non -homogeneous Eq. (63). By neglecting 
e 2 d we have for a particular solution of Eq. (63) 


ot 


e 


€ 

(I 7 * 7 !) 


[c 


sinT + d cos t] - 


e 2 bc sin 2t 
2(a 2 


(72) 


In the last expression we can see that a 2 = 1 corresponds to a resonance. 
Our results constitute an extension of Beletskies study of small eccentric 
oscillations of a satellite under pure gravity torque (Ref. 6, p. 41). The 
general solution of Eq. (63) is then the sum of the general solution of the 
Mathieu equation and the particular integral, Eq. (72). 


CONCLUSION 

In this paper we have presented an analytical study of the longitudinal dyna- 
mics of a thrusting, lifting, orbital vehicle in a nearly circular orbit. Ex- 
plicit expressions for the elements of the orbit were derived and the behaviors 
of the variations of these elements were correctly predicted. It was shown 
that for large perturbations the second order gradient effect of the air mass 
density must be included. Explicit expressions for the period and damping 
of the angle-of-attack mode were derived. It was shown that a resonance 
effect was not present for a circular orbit. A resonance effect was displayed 
by a study of the forced eccentricity oscillations and the critical altitude for 
resonance was obtained by solving a very simple equation. The analytical ex- 
pressions are in excellent agreement with an independent numerical analysis 
at all altitudes. 
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APPENDIX A 


Characteristics of the Atmosphere 

All numerical computations are based on the U.S. Standard 
Atmosphere, 1962*. However, for ease of computation the approximate 
inverse polynomial representation of this atmosphere that appears in the 
U.S. Standard Atmosphere Supplements, 1 96 6 1, is assumed to be exact. 

The inverse polynomial is of the form 

p/p o = l/[A 0 + A 1 Z+- • •+A 11 Z 11 ] 4 (A- 1) 

s 

where Z is the geometric altitude above the standard geoid (637 8. 17 km 
radius) in kilometers. The coefficients A^ are given in Table A-l. 

The polynomial approximation is valid for the altitude range 0-200 
km. Compared to the Standard Atmosphere, the approximation differs by 
less than 5% (see Fig. A-l) in this altitude range. 

Error % * 

0 40 eo (20 260 200 

Fig. A-l Error in Representing M 62 M 
Standard Atmos. By Poly- 
nomial Representation, f 

The non-dimensional density gradients defined in Eq. (10) of the 
text are plotted in Fig. A- 2. For comparison the first density gradient 

* U.S. Standard Atmosphere, 1962, prepared under sponsorship of ESSA, 
NASA, and USAF. 

|U.S. Standard Atmosphere Supplemen ts, 1966 , prepared under sponsor- 
ship of NASA, USAF and U.S. Weather Bureau. 

Both publications are available from the Superintendent of Documents, 

U.S. Government Printing Office, Washington, D.C. 20402. 

f lbid . , p. 68. 
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(o-j) obtained by numerical differentiation of the tabulated Standard 
Atmosphere data is included in the figure. 

Table A-l 

COEFFICIENTS FOR DENSITY POLYNOMIAL 
Altitude range 0-200 lan 
= 1.2250 Kilograms/m 3 A_. [=] W 


3 

A. 

3 

0 

0. 1000000000 

E 01 

1 

0. 3393495800 

E-01 

2 

-0.3433553057 

E-02 

3 

0. 5497466428 

E-03 

4 

-0.3228358326 

E-04 

5 

0. 1106617734 

E-05 

6 

-0. 2291755793 

E- 07 

7 

0.2902146443 

E-09 

8 

-0.2230070938 

E-ll 

9 

0. 1010575266 

E- 13 

10 

-0. 2482089627 

E- 16 

11 

0. 2548769715 

E- 19 
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APPENDIX B 


Characteristic of the Vehicle 


The characteristics of the vehicle geometry and the aerodynamic 
stability derivatives used in all numerical calculations are those used by 
Etkin 1 and Rangi 4 and are typical for a slender body, a cone or wedge of 
3° semiangle. The values were derived from the simple Newtonian im- 
pact theory for moderate angles of attack. They are tabulated below. 

Table B-l 

GEOMETRIC VEHICLE PARAMETERS 
k = 6 ft. L = 50 ft. 

y 

k 0 = -0.94 W/S = 30 psf . (sea level) 


Table B-2 


AERODYNAMIC STABILITY DERIVATIVES 

8C 


C T = 0.05 


D 0 


= 0.0133 


m. 


m 




= -0. 028 


a 


dC 

— ! - = 0.329 C^ 
da D 


a 


dC 


D 


Q =0.15 C 
da m 


9C 


m 


da 


= -0. 0548 
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APPENDIX C 


Effects of Second Order Atmospheric Mass Density Gradient 

In pages 13 and 14 we have discussed the need for inclusion of the 
second order atmospheric mass density gradient (terms containing c r 2 ) 
for an accurate analysis. This effect can be dramatically illustrated by 
a numerical analysis. 

Figures C-l and C-2 respectively represent the variations of the 
radial distance and the flight path angle as the time varies. They are 
reproductions of computer generated plots with different atmospheric 
mass density laws. The set is generated at an altitude of 300, 000 ft, 
where the density gradients are the targets, with an initial non-dimen- 

— 3 . — 3 

sional speed decrease of 10 (e = Au - -10 ). This corresponds to a 

perturbation of 25 ft/sec. The equations used are the uncoupled phugoid 
equations (Eqs. 7, p. 5). The solid lines represent the time histories of 
the elements of the orbit obtained by a numerical integration of the un- 
coupled non-linear equations, Eqs. 7. One curve is generated with the 
exact density law, that is the 44th degree inverse polynomial representa- 
tion of the 62 Standard Atmosphere. Another curve is obtained by keep- 
ing only the second order gradient term. The dotted line represents the 
linear analytical solutions of the uncoupled motion (Eq. 18) with the value 
for a to the order . For comparison the numerical solutions using 
Etkin T s linearized and coupled equations are also plotted. It is clear from 
the graphs that our linear analytical solutions and Etkin T s linear numerical 
solutions are nearly identical. But they do not compare well after 1 / 4 of 
a revolution with the exact solution using the exact atmospheric mass 
density law. The result can be improved by including higher order den- 
sity gradients. This is obtained by using our second order solution 
(Eq. 22). The discrepancies are much less at higher altitude where 
atmospheric drag is small and at lower altitude where the period is short 
and practical perturbations are small. The time scale taken for the plots 



Non Dimensional Radius 



Fig. C-l Variations of Radial Distance 







is about 2 linear periods. Referring to Fig. C-2 the flight path angle 
time history shows that the exact integration curve displays less spiral 
mode and less damping than the linear solution and further shows that the 
exact phugoid oscillations have smaller period. In fact the first cycle 
takes 1.95 X 10 3 secs versus the linear phugoid period of 2. 51 X 10 3 secs. 
The second and third cycles take even less time and show that the exact 
motion will complete 3 cycles for the two linear periods. 

A complete numerical analysis shows that below 100, 000 ft or above 
400, 000 ft linear solutions are accurate. In between there is a definite 
requirement for the inclusion of higher order atmospheric mass density 
gradients. 
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